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Time  Lapse  Simulation  of 
Interrelated  Weather  Conditions 


1.  INTRODUCTION 

Climatic  records  give  the  frequencies  of  weather  elements  such  as  ceiling  and 
visibility,  showing  the  cyclic  variations  with  time  of  day  and  season  of  the  year. 
Weather  records  of  hourly  observations  provide  sequences  of  events  or  of  temporal 
changes  of  ceiling  and  visibility.  However,  there  is  a  need  for  simulation  of  these 
events  as  a  stochastic  process  that  takes  into  account  such  attributes  of  the  weather 
as  the  association  of  the  visibility  with  the  ceiling  and  the  degree  of  persistence  of 
both  elements  from  hour  to  hour. 

Figure  1  is  a  graphical  record  of  actual  Boston  ceiling  and  visibility  from 
midnight  (Z),  22  Dec  1982  to  midnight  (Z)  48  h  later.  The  first  day  was  cloudless 
and  clear,  with  both  ceiling  and  visibility  "unlimited."  On  the  second  day  clouds 
moved  in,  the  ceiling  lowered,  and  visibility  decreased. 

Figure  2  shows  the  same  information  for  the  celling  on  normal  probability 
paper.  There  are  several  advantages  in  plotting  the  sequence  of  cloud  ceilings  on 
Figure  2  rather  than  Figure  1.  The  dashed  lines  in  Figure  2  show  the  climatic 
cumulative  frequencies  throughout  the  day  of  ceilings  from  500  ft  to  above  30,  000  ft. 
This  graph  provides  the  additional  information  that  the  low  ceilings  of  500  ft  or  less 
have  an  a  priori  probability  of  2  to  3  percent  in  December,  which  is  roughly  for 
20  h  during  the  month.  Likewise,  Figure  3  shows  the  climatic  frequencies  for 
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Figure  1.  The  Hourly  Sequence  of  Ceiling  Height  and  Visi 
bility  at  Boston,  Mass.,  From  1900  EST,  21  Dec  1982  to 
1900  EST,  23  Dec  1982 
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Figure  2.  The  Hourly  Sequence  of  Ceiling  Height  at  Boston, 
Mass.,  From  1900  EST.  21  Dec  1982  to  1900  EST,  23  Dec 
1982.  Dashed  lines  give  the  climatic  frequencies  of  ceiling 
as  a  function  of  time  of  day 
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Figure  3.  The  Hourly  Sequence  of  Visibility  at  Boston, 
Mass.,  From  1900  EST.  21  Dec  1982  to  1900  EST,  23  Dec 
1982.  Thin  dashed  lines  give  the  climatic  frequencies  of 
visibility  as  a  function  of  time  of  day 


visibility,  in  comparison  with  the  actual  sequence  during  the  2 -day  period.  As 
may  be  seen  in  subsequent  examples  in  this  paper,  an  improvement  or  a  deteriora¬ 
tion  in  the  .  eiling  and  visibility  could  be  a  function  of  the  time  of  day. 

The  kind  of  record  just  described  is  not  readily  available.  Even  if  it  were, 
using  it  in  tests  to  provide  a  realistic  sequence  of  events  would  still  be  difficult. 
During  a  period  as  long  as  50  years,  many  meteorological  situations  arise,  yet 
they  do  not  exhaust  all  the  possibilities. 

This  report  presents  a  stochastic  procedure  for  simulating  the  changes  over 
time  in  ceiling  and  visibility  that  are  characteristic  of  a  real  climate.  It  is 
assumed  that  the  changes  in  ceiling  and  visibility  are  a  Markov  process  (as  de¬ 
scribed  in  Section  3).  There  is  a  significant  correlation  between  present  ceiling 
and  present  visibility.  When  there  is  a  time  lag  between  the  observations  of 
ceiling  and  visibility,  the  correlation  coefficient  is  reduced.  The  effect  is  explored 
in  this  report. 
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2.  EQUIVALENT  NORMAL  DEVIATES  (END) 


Many  weather  elements,  such  as  ceiling  and  visibility,  have  distributions  that 
are  difficult  to  treat  statistically.  Ceiling  frequencies  are  usually  given  in  terms 
of  categories  of  height  above  the  ground,  and  there  is  usually  a  substantial  prob¬ 
ability  of  occurrence  of  "unlimited"  ceiling;  that  is,  clear  sky  or  scattered  clouds 
Modeling,  however,  has  been  reasonably  successful.  Bean  et  al*  have  used  Burr 
curves  for  the  cumulative  distribution  of  ceiling  heights  in  the  form 

FW-l-  [i.  (i)3] 

where  x  is  ceiling  height,  in  feet 
ceiling  equal  to,  or  less  than,  x, 
to  provide  the  best  least  squares 
January,  from  12  to  14  LST,  the 

a  =  1.  16779 
b  =  0. 192682 
c  =  1000  . 

which  give  the  solid  curve  plotted  in  Figure  4.  The  RUSSWO*  informat'  n  on  the 

cumulative  frequencies  of  ceiling  height  are  shown  (by  x's)  on  Figure  4.  The 

model  fits  well,  but  requires  some  caution  in  application,  because  unless  it  is  the 

right  tool  for  the  problem,  the  answer  it  gives  might  be  misleading.  For  example, 

the  ceiling  is  unlimited  with  an  observed  frequency  of  0.  445.  This  is  comparable 

to  the  formula  estimate  of  probability  of  ceiling  above  32,  000  ft.  Ceiling  less  than 

100  ft  should  be  categorized  with  a  probability  of  0.  005. 

Visibility  frequencies,  likewise,  have  been  fitted  by  an  idealized  model. 

2 

Somerville  et  al  use  the  Weibull  distribution  to  give  the  cumulative  probability 
of  visibility  (x)  as 

F(x)  =  1  -  exp  (-ax^)  (2) 


a,  b.  c>  0  (1) 

or  meters,  F(x)  is  the  cumulative  probability  of 
and  a,  b,  and  c  are  model  parameters  estimated 
fit.  For  Bedford,  Mass.,  for  example,  in 
values,  as  given  by  Bean,  are 


% 

RUSSWO  stands  for  "Revised  Uniform  Summaries  of  Surface  Weather  Obser¬ 
vations,  "  published  by  USAF  Environmental  Technical  Application  Center  for 
several  hundred  stations  around  the  world. 


1.  Bean,  S.J.,  Somerville,  P,  N. ,  and  Heuser,  M.  (1979)  Some  Models  for 

Ceiling,  Scientific  Report  No.  7,  Contract  F19628-77 -C-0080,  ADA078033. 

2.  Somerville,  P.  N. ,  Bean,  S.J.,  and  Falls,  F.  (1979)  Some  Models  for  Vis¬ 

ibility.  Scientific  Report  No.  3,  Contract  F19628-77-C-0080,  ADA075490. 


12 


UNLIMITED 


I  10  I02  I03  I04  I05 

CEILING  HEIGHT-(ft) 


Figure  4.  The  Cumulative  Climatic  Frequency 
of  Ceiling  Height  at  Bedford,  Mass. ,  in  the 
Month  of  January,  12-14  EST.  The  X's  mark 
the  frequencies  as  given  in  RUSSWO  tables. 

The  solid  line  is  the  "best"  model  fit  by  a  Burr 
curve 


where  a  and  p  are  parameters.  For  Bedford,  Mass.,  in  January,  from  12  to  14 
LST,  the  values  are: 

o  =  0.06906 

P  =  0.8186  , 

which  give  the  curve  as  plotted  in  Figure  5.  The  RUSSWO  data  again  reveal  that 
the  model  fits  well,  but  visibility  less  than  1/4  mile  should  be  categorized  with  a 
probability  of  2  percent.  Visibility  greater  than  10  miles  should  be  categorized 
with  a  probability  of  62  percent. 

Any  probability  F(x)  as  given  in  Eq.  (1)  or  (2),  corresponds  to  an  equivalent 
normal  deviate  (END),  symbolized  by  y(0,  1).  The  latter  is  a  variable  with 
Gaussian  distribution,  mean  value  of  0,  and  standard  deviation  1.  0.  Symbolically, 


•T'T 


Figure  5.  The  Cumulative  Climatic  Fre¬ 
quency  of  Visibility  at  Bedford,  Mass. ,  in 
the  Month  of  January,  12-14  EST.  The  X’s 
mark  the  RUSSWO  frequencies.  The  solid 
curve  is  the  "best"  model  fit  by  a  Weibull 
distribution 


y  2 

F(x)  =  — ^  f  exp  -  -X-  dy  ,  (3) 

-oo 

which  thus  defines  the  END  (y).  The  latter  is  implicitly  found  for  each  x. 

3 

To  a  highly  satisfactory  approximation,  (error  <  0.  003) 


y  =  k 


t 


ao  +  alt 
1  +  bxt  +  b2t2 


where 


a  =  2.30753 
o 

a:  =  0.27061 

3.  NBS  (1964)  Handbook  of  Mathematical  Functions.  Applied  Mathematics  Series, 
No.  55,  Government  Printing  Office,  Washington,  D.  C. 


bl  =  0.99229 
b2  =  0.04481 


and 


In  -4j  when  p  =  F(x)  <1/2 
P 

(4) 


In - - — 5-  when  p  =  F(x)  >  1/2  . 

(1  -  p^) 

A  scale  of  END,  (y)  appears  alongside  the  scale  of  F(x)  in  both  Figures  4  and  5. 

In  this  report,  ceiling  height  and  visibility  are  modeled  in  terms  of  ENDs 
because  they  have  varying  averages  and  medians  throughout  the  day,  and  will  have 
generally  uneven  distributions.  But  the  END  of  the  ceiling  or  visibility  will  have 
a  symmetrical  Gaussian  distribution  with  a  constant  mean  or  median  of  zero.  In 
Figure  3  it  can  be  seen  how,  for  example,  visibility  less  than  5  miles  has  varying 
probability  throughout  the  day.  The  END  varies  directly  with  this  probability. 
Correlations  in  this  paper  are  found  between  the  END'S  of  ceilings  at  differing 
times  of  day,  or  between  the  END'S  of  ceiling  and  visibility,  with  or  without  time 
lag. 


k  =  1,  t 

or 

k  =  -1, 


3.  THE  ORNSTEIN  UHLENBECK  (0-U)  PROCESS 

A  Markov  process  is  defined  as  "a  stochastic  process  such  that  the  conditional 

probability  distribution  for  the  state  at  any  future  instant,  given  the  present  state, 

4 

is  unaffected  by  any  additional  knowledge  of  the  past  history  of  the  system". 

The  Ornstein-Uhlenbeck  (O-U)  process  is  one  kind  of  Markov  process  in  which 
the  value  at  a  future  instant  t  +  6t,  of  a  normally  distributed  variable  y  is  linearly 
related  to  the  value  at  the  present  instant  t,  and  the  correlation  coefficient  between 
present  and  future  values  decays  exponentially  with  the  time  interval  it  between 
them.  Mathematically,  this  may  be  stated  as  follows: 


4.  Kendall,  M.G.,  and  Buckland,  W.R.  (1971)  A  Dictionary  of  Statistical  Terms, 
Hafner  Publishing  Co. ,  New  York,  .  ~ 
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where  the  correlation  coefficient  p  is  given  by 


p  =  exp  (-6t/  r  )  (6) 

where  t  is  a  parameter  with  the  dimension  of  time,  named  the  Relaxation  Time. 

5  6 

by  Keilson  and  Ross,  and  described  by  Gringorten.  The  random  number  p  is 
normally  distributed  and  is  the  random  contribution  to  the  process  of  change. 

A  time  sequence  of  y's  at  regular  intervals  (6t),  generated  by  the  model  of 
Eq.  (5),  will  simulate  the  conditions  of  a  weather  element  in  terms  of  its  END. 

A  step-by-step  program  to  accomplish  this  is  given  in  Appendix  A. 

Figure  6a  was  drawn  for  a  relaxation  time  t  of  20  h,  which  corresponds  to  a 
realistic  hour-to-hour  correlation  coefficient  of  0.  95.  Figure  6a  presents  the 
simulation  of  sky  cover  for  a  partly  cloudy  to  overcast  48-h  period.  The  dashed 
lines  show  climatic  frequencies  of  clear,  scattered,  broken  and  overcast  through¬ 
out  a  day  in  August  at  Bedford,  Mass. ,  deliberately  chosen  because  of  a  large 
diurnal  effect.  Figure  6b  illustrates  a  sequence  in  which  serial  correlation  is 
weaker  (0.82)  with  T  =  5  h.  Still,  with  this  much  persistence  an  overcast  might 
remain,  with  few  breaks,  for  20  h.  There  is  a  clearing  at  the  most  likely  time  of 
the  24-h  period.  Figure  6c  illustrates  a  Sequence  in  which  serial  correlation  is 
stronger  (0.98)  with  T  -  50  h.  Finally,  Figure  6d  illustrates  a  sequence  in  which 
hour-to-hour  correlation  is  reduced  virtually  to  zero,  with  r  =  0.  1  h,  producing 
rapid  changes  and  short  periods  of  all  sky -cover  conditions. 


5.  Keilson,  J. ,  and  Ross,  H.  F.  (1979)  Gaussian  Markov  Related  Variates  for 

Meteorological  Planning.  Final  Report.  Contract  F19628-78-C-0158 
AD  A081382.  - 

6.  Gringorten,  1. 1.  (1982)  The  Keilson-Ross  Procedure  for  Estimating  Climatic 
Probabilities  of  Duration  of  Weather  Conditions.  AFGL-TR-82-0116 
AD  A1 19860. 


Figure  6a.  Simulation  of  a  48-h  Sequence 
of  END's  in  an  O-U  Stochastic  Process, 
With  a  20 -h  Relaxation  Time.  The  dashed 
lines  show  climatic  frequencies  of  clear, 
scattered,  broken,  and  overcast  through¬ 
out  a  day  in  August  at  Bedford.  Mass. 


4.  TWO  END’S  IN  A  JOINT  MARKOV  PROCESS 


The  above  procedure  gives  a  time  lapse  simulation  of  a  single  weather  ele¬ 
ment.  The  next  task  is  to  provide  the  algorithms  for  time  lapse  simulation  of  two 
elements  that  are  interrelated.  Two  such  weather  conditions  (X.,  X..)  may  be  the 
two  sky  covers  at  two  neighboring  stations,  or  they  may  be  ceiling  and  visibility 
at  one  station.  Whatever  they  be,  in  this  paper  they  are  represented  by  their 
END'S  (y.,  v..),  each  of  whose  means  is  zero,  standard  deviation  1.0,  and  whose 
probability  distribution  is  Gaussian. 

The  two  END'S,  y . (t)  and  yTt)  corresponding  to  the  two  events  (X.,  X.)  at  time 
(t)  are  subject  to  change  in  the  time  interval  <6 1)  in  a  Markov  process,  given  as 
follows: 


y . (t  +  6t)  =  a.j  •  y.(t)  +  a^  •  .Vj(t)  +  b.  •  p.(t  +  6t) 
y,(t  +  6t)  =  a.j  •  y.(t)  +  a..  •  y^t)  +  b  .  p.(t  +  6t) 


(7) 


where  rjj.  hj  ace  normally  distributed  and  random,  except  for  their  interrelation¬ 
ship.  The  a's  are  partial  regression  coefficients  and  the  b's  are  of  such  magni¬ 
tude  that  the  normality  of  y.,  y^  is  preserved.  The  a's  and  b's  need  to  be  deter¬ 
mined  in  terms  of  correlation  coefficients,  which  are  derivable  from  the 
climatology  of  the  station  or  stations. 

Because  of  the  normality  of  the  END'S  the  correlation  coefficients,  by 
definition,  are: 


rij =  Eiyi<t>  *  yj^l  =  EIyi^  +  5t*  *  yi<1  +  6tM 

where  the  symbol  E[  ]  denotes  the  expected  value  of  the  quantity  in  brackets,  and 
r.j  is  the  correlation  coefficient  between  the  elements  without  time  lag.  For  time 


lag  (6t): 

PU  » 

Efyjit)  • 

yi(t 

+  St)] 

Pij  = 

E  [y.(t)  • 

yj(t 

+  «5t)] 

Pii  = 

E[y.(t)  • 

yl<‘ 

+  «5t)] 

Pjj 

E[y.(t)  • 

yj(t 

+  5t)] 

(9) 
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Also,  the  variances  are: 


E  yf(t)  =  E  y?(t)  =  E  nf  =  E  n*  =  1  . 

The  mean  and  variance  of  each  END  do  not  vary  with  time. 

From  Eqs.  (7)  and  (8),  after  squaring  both  sides  of  each  equation: 

i  =  4  +  afj  +  2aiiaijrij  +  b? 


1  =  a?.  +  a2  +  a . .  a .  ■  r . .  +  b2 
33  31  33  13  3 


The  derivation  of  these  equations  takes  advantage  of  the  fact  that  17.  and 
obtained  independently  of  y^  and  y^. 

Again,  from  Eqs.  (7)  and  (9): 


f>ii 

=  au 

+  ay 

r. . 

*3 

p  *  * 

'  aji 

+  a. . 

r. . 

33 

l3 

P3i 

=  aii 

r..  + 

a;. 

l) 

U 

P  ■  • 

=  a. . 

r\.  + 

a. . 

j3 

Ji 

i) 

33 

Solving 

for  th 

e  a's: 

aii 

=  (Pj, 

*  pji 

r..)/(l 

2  , 

a.. 

=  (P  .  . 

*  P;; 

r..)/(l 

-r2) 

i) 

Jl 

^  ll 

13 

l.l 

a)i 

=  <Pii 

'  P33 

rV/a 

- 

a . . 
.1.1 

=  (P.ij 

■p*i 

r..)/(l 

t.l 

-  r2l 
1) 

3  are 


From  Kq.  (10): 


There  remains  the  question  of  the  interrelation  of  the  stochastically  produced 
values,  rjjd  +  fit)  and  ^ (t  +  fit).  Symbolizing  the  correlation  coefficient  between 
them  as  h.  we  find  from  Eq.  (7): 


b.bj  .  h  =  E  {yt(t  +  fit)  -  a..  •  yt<t)  -  a^  •  y^t)}  • 
{yj(t  +  fit)  -  a...  •  y.(t)  -  a...  .  y^t)} 


Hence, 


h  =  {r..(l  -  a.,  a..  -  a.  a..)  -  (a.,  a.j  +  a{j  a.^J/b.  b. 


(14) 


If  two  random,  normal  numbers  r)j(t  +  fit),  rjj <t  +  fit)  are  selected,  then: 


r).(t  +  fit)  =  h  •  r).(t  +  fit)  +  V  1  -  h2  •  nj(t  +  fit) 


(15) 


Thus,  all  the  values  are  obtainable  for  the  solution  of  Eq.  (7)  if  the  correlation 
coefficients  are  provided.  By  substitution  of  y^(t  +  fit),  y^(t  +  fit)  for  yj(t),  (t) 
the  equations  can  be  solved  for  the  next  pair  of  yt  and  y^.  By  such  iteration,  the 
simulation  of  a  joint  sequence  of  values  of  yj,  y^  is  obtainable.  The  computer 
programming  of  this  process  is  given  in  Appendix  B. 

7 

A  variation  of  the  above  solution  was  obtained  by  Maj.  R.  C.  Whiton  at 
USAF/ETAC,  Scott  AFB,  Ill.  Whereas  Eq.  (7)  presupposes  that  the  later  value 
of  each  of  the  ceiling  and  visibility  is  dependent  on  both  the  previous  ceiling  and 
the  previous  visibility,  Whiton's  equation  for  the  later  ceiling  assumed  its 
dependence  on  the  current  ceiling,  but  not  on  the  current  visibility;  likewise  for 
the  later  visibility.  Thus  the  process  of  change  in  either  ceiling  or  visibility 
satisfies  the  definition  of  an  O-U  process.  This  is  equivalent  to  setting  the  lag 
correlation  equal  to  the  product: 


p  . .  =  r. .  •  p  . . 


r . .  •  o  ■■ 
ij 


(16) 


Whiton's  simplification  will  often  prove  effective.  However,  for  generality  it 
might  be  better  to  avoid  this  assumption  and  let  the  climatic  data  provide  numbers 
for  p.j,  p^.  Some  sources  of  the  correlation  coefficients  are  given  and  discussed 
in  the  following  sections. 


7.  Whiton,  R.C.  (1982)  Environmental  Simulation  Modeling,  AWS/ETAC, 
Scott  AFB,  Ill. 
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5.  ESTIMATING  CORRELATION  COEFFICIENTS  BETWEEN  CEILING  AND  VISIBILITY 


A  measure  of  the  intercorrelation,  without  time  lag,  between  ceiling  and 
visibility  is  obtainable  from  a  RUSSWO.  Table  1,  for  example,  presents  the 
joint  probabilities  (Pfe  C,  >V)),  as  interpolated  from  the  RUSSWO  tables  for 
Bedford,  MA,  midnight,  in  January,  April,  July,  and  October.  The  probabili¬ 
ties,  as  selected,  are  the  same,  in  pairs,  for  ceiling  and  for  visibility.  The 
corresponding  END's  are  also  shown  in  Table  1.  "Tables  of  the  Bivariate  Normal 
Distribution"  (NBS,  1957)  make  it  possible  to  find,  by  interpolation,  the  inter¬ 
correlation  coefficient  (p),  given  the  probabilities  Pfe  C),  P(z  V),  or  their 
END's,  together  with  their  joint  probability  P(s  C,  a  V).  Correlation  that  is 

9 

estimated  from  these  probabilities  is  known  as  tetrachoric  correlation. 


Table  1.  Estimates  of  the  Joint  Cumulative  Probability  of  Ceiling,  F(s  C), 
and  Visibility,  P(>V),  at  Bedford,  Mass.,  at  Midnight.  The  conditions  of 
ceiling  and  visibility  are  those  corresponding  to  their  percentiles:  60,  70, 
80,  and  90  percent 


END 

Month 

January 

April 

July 

October 

mam 

MSB 

wbsm 

0.48 

0.43 

0.42 

■SI g: 

I 

0.  60 

0.  56 

0.  54 

0.8 

0.  84 

I 

0.74 

0.70 

0.73 

0.  9 

1.28 

MSM 

0.  86 

0.85 

0.  86 

A  computerized  solution  of  the  intercorrelation  coefficient  may  be  found  to 
satisfy  the  relation  upon  which  the  NBS  tables  are  based: 


P(>x,  >y)  = 


2*yji 


// 


g2  +  n2  -  2 png 

2d  -  p2) 


d?  •  drj 


(17) 


Programming  such  an  equation  on  a  desk-top  computer  is  difficult  and  provides 
only  approximate  solutions.  The  computer  operation  is  slow  because  it  requires 
trial -and-error  iterations. 

8.  NBS  (1957)  Tables  of  the  Bivariate  Normal  Distribution  Function,  Applied 

Mathematics  Series.  No.  50,  Government  Printing  Office,  Washington,  D.G. 

9.  Brooks,  C.E.P.,  and  Carruthers,  N.  (1953)  Handbook  of  Statistical  Methods 

in  Meteorology,  Her  Majesty's  Stationery  Office,  London. 


A  reasonable  estimate  of  the  tetrachoric  correlation  coefficient  has  been 


presented,  as  follows: 


xy 


sin 


*_  -s/ad"  -%/bF 
2 


%/ad  +  %/bc" 


(1 


where 


a  =  P(>  X,  >  Y) 
b  =  P(>  X)  -  a 
c  =  Pfe  Y)  -  a 
d  =  1-a-b-c 


The  results  of  using  this  formula  are  shown  in  the  parentheses  (Table  2),  for 
easy  comparison  with  r.ne  results  of  applying  a  computerized  approximation  of 
the  bivariate  normal  distribution.  Since  there  is  little  or  no  significant  differ¬ 
ence  in  the  results,  Eq.  (18)  is  favored  because  it  is  a  simpler  method  for  find 


ing  the  correlation  coefficient  (r  ). 

xy 


Table  2.  Bivariate  Normal  Distribution  Estimates  of  (rpy)>  the  Tetrachoric 
Correlation  (Zero  Time  Lag)  Between  Ceiling  and  Visibility,  at  Several  END 
Levels,  for  Bedford,  Mass.  Midnight  Values.  Figures  in  parentheses  were 
obtained  by  using  the  model  of  Eq.  (18),  from  Brooks  and  Carruthers" 


P(2r  C) 

P(>V) 

END 

January 

Month 

April 

July 

October 

0.  6 

0.25 

0.  66 

0.72 

0.47 

0.40 

(0.  67) 

(0.72) 

(0.45) 

(0.39) 

0.7 

0.  53 

0.78 

0.  78 

0.  50 

0.40 

(0.76) 

(0.76) 

(0.54) 

(0.40) 

0.8 

0.84 

0.  91 

0.84 

0.  59 

0.84 

(0.  92) 

(0.88) 

(0.  65) 

(0.84) 

0.  9 

1.28 

0.  93 

0.90 

0.  80 

0.85 

(0.  98) 

(0.  89) 

(0.82) 

(0.89) 

Using  the  formula  of  Eq.  (18),  the  tetrachoric  correlation  coefficient  was 
found  for  the  RU3SWO  frequencies  of  ceiling  (C)  and  visibility  (V)  and  joint  fre¬ 
quencies  at  Bedford,  Mass. ,  in  January,  at  midnight  and  at  12  noon  (Table  3). 
There  are  no  entries  for  ceiling  less  than  2,000  ft  or  for  visibility  less  than  2 
miles  because  the  samples  were  too  small  to  give  meaningful  figures. 


Table  3.  Non-Lag  Tetrachoric  Correlation  Coefficients  Between  the  END's  of 
Ceiling  (C)  and  Visibility  (V)  at  Bedford,  Mass.  ,  in  January,  Derived  From 
RUSSWO  Tables.  In  each  box  the  upper  figure  is  the  midnight  value,  the  lower 
figure  the  noontime  value 


V 

>  10 

>  6 

>  4 

>  2 

c 

P(>C) 

PteV) 

0.  592 

0.  623 

0.735 

0.746 

0.831 

0.818 

0.887 

0.884 

>  30,000 

0.  507 

0.64 

0.72 

0.  81 

0.82 

0.445 

0.78 

0.  89 

0.94 

1.00 

=>  10,000 

0.  597 

0.70 

0.76 

0.84 

0.84 

0.  599 

0.  80 

0.87 

0.  94 

1.00 

>  5,000 

0.  684 

0.75 

0.  82 

0.86 

0.  88 

0.675 

0.80 

0.  86 

0.  92 

0.  98 

2  2,000 

0.805 

0.  92 

0.  91 

0.92 

0.  94 

0.805 

0.95 

0.  93 

0.  94 

0.  97 

Unfortunately,  the  tetrachoric  correlation  coefficient  shows  considerable 

dependence  on  the  specific  values  of  ceiling  and  visibility,  or  alternatively  on 
the  probabilities  of  the  events.  It  can  be  expected  that  low  ceilings  and  visibili 
ties  will  be  highly  correlated,  both  having  low  probabilities  of  occurrence.  We 
next  examine  this  effect  further. 


6.  TWO  JOINT  WEATHER  SEQUENCES 

By  the  procedure  of  Section  4,  we  are  able  to  generate  pairs  of  values  of 
END's,  for  example,  one  for  ceiling  and  one  for  visibility.  For  practical  pur¬ 
poses  these  should  be  transformed  into  units  of  ceiling  and  visibility  at  one  sta¬ 
tion,  or  into  the  units  of  ceiling  at  two  neighboring  stations. 

O 

An  END  (y)  can  be  transformed  into  its  probability,  which  is  also  the 
probability  of  the  original  variable  (X).  Thus,  using  the  notation  P(£  X)  for  the 
cumulative  probability  of  X,  we  are  given  (with  3 -decimal  accuracy): 


P(SX)  =  P(<y)«i+m  [2(1  +  c  |y|  +  c2y‘ 


,2  +  c3 1  y3  | 


where 


Cj  “  0.  196854 

c2  *  0.  115194 

c3  =  0.  000344 

c.  =  0.019527 
4 


t  -  1,  m  =  -1  for  y  >  0 
f  =  0,  m  =  1  for  y  <  0 

If  the  Burr  curve  *  is  accepted  as  the  model  for  the  distribution  of  ceiling  height 
(C)  then  C  is  given  by 


C 


c[(l  -  p)‘l/b  - 


(2-0) 


where  p  =  P(i  C);  a,  b,  and  c  are  parameters.  Bean  et  al*  give  tables  of  values 

for  these  three  parameters  to  yield  answers  for  ceiling  heights  in  feet. 
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If  the  Weibull  distribution  is  accepted  as  the  model  for  visibility  (V)  then 
V  is  given  by 


V  =  fn  (1  -  P)j  (21) 

2 

where  p  «  Pfe  V),  and  a  and  |3  are  parameters.  Somerville  et  al  give  tables  of 
values  for  these  two  parameters  to  yield  visibility  in  miles. 

The  procedure  for  generating  a  stochastic  sequence  of  joint  events  is  linked 
to  the  previously  outlined  procedure  (Section  4)  in  the  following  order: 

At  each  step  in  the  sequence,  pairs  of  random  normal  numbers  are  selected 
and  used  with  the  known  correlation  coefficients  for  the  determination  of  the 
END's  (Eq.  (7)].  Each  END  (y)  is  then  transformed  into  the  corresponding 
probability  p  =  Pfe  y)  by  Eq.  (19).  Then  each  p  is  transformed  into  the  weather 
element  itself,  into  ceiling  height  by  Eq.  (20),  or  into  visibility  by  Eq.  (21). 

The  computer  program  is  in  Appendix  C. 


7.  MULTIPLE  JOINT  EVENTS  IN  A  MARKOV  PROCESS 


The  above  derivation  could  be  generalized  for  n  weather  elements  in  the 
form: 


^  Xj  b^  •  ni  for  i  =  l.n  (22) 

i 

where  the  y's  are  the  predictands  and  the  x's  are  the  predictors.  Where  A  is  a 

T  1 

column  vector,  A.  is  a  row  vector  of  n  coefficients: 


1i2 . ain) 


(23) 


Squaring  Eq.  (22)  and  taking  expected  values: 


b.  =  ^ 1  -  A?  C  A  for  i  =  l.n  (24) 

where  C  is  the  matrix  of  correlation  coefficients  (r.J  between  the  n  predictors 
without  lag: 


!*  r12 . rln 

r2  1 


Eq. 


(22)  represents  n  equations,  each  of  the  form 
ailXl  +  -"  +  ain  xn  +  bi  r,i  =  ?i  f0ri=1'n 


(25) 


(26) 


To  solve  for  the  n  coefficients  (a.j,  ....  a.n>  we  find  n  equations  by  multiplying 
Eq.  (26)  successively  by  Xj,  ....  and  then  taking  expected  values.  This  pro¬ 
cedure  introduces  the  correlation  coefficients  (r.^)  between  the  predictors  without 
time  lag,  and  the  correlation  coefficients  (p^)  between  the  predictand  (yj  and 
each  of  the  predictors  (x.,  j  =  l.n).  Thus  the  n  equations  become 


From  Eq.  (22): 


b.  n.  =  v-  -  a.„x,  ...  -a.  x 
1  '1  Ji  ll  1  in  n 


(32) 


b.  n.  =  y.  -  a.,x,  .  . .  -a.  x 
j  ']  il  1  )n  n 


If  hy  is  the  correlation  coefficient  between  and  then 


b. 


ibjhij  -  rij  -  E  aikpjk  *  E  ajk  pik+  E  a 

k=l  k=l  k,  l 


ik  '  aji  rkf 


(33) 


which  defines  the  (  2  )  correlation  coefficients  between  the  17's.  [(  2  )  represents 
the  number  of  possible  combinations  of  n  things,  taken  two  at  a  time. )  To  find 


values  for  q.,  (j  =  l,n),  for  use  in  Eqs.  (19),  first  choose  n  random,  independent 
numbers 


($..  i  =  1,  n)  . 


The  relations  of  the  q's  to  the  ? 's  are  of  the  form: 


°11?1  +  +  *ln?n  = 


"21  ?1  +  •••  +  °2n?n  =  *2 


(34) 


+  a  £  =  r) 

nn  n  n 
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The  task,  now,  is  to  find  suitable  values  for  the  a's.  There  are  n  values  for 

2 

the  a's  but  only  n  +  (  2  )  e9uati°ns  are  obtainable  by  setting  Eq.  =  1,  and  by 
finding  the  (  £  )  values  of  h.^  from  Eq.  (33).  A  certain  arbitrariness,  therefore, 
is  preserved  in  assigning  values  to  some  of  the  a's. 

For  n  =  2,  or  for  only  two  weather  elements,  as  seen  previously,  if  we  can 
set  on  =  1,  it  follows: 


belief  that  modeling  the  probability  of  fractional  cover  of  areas  of  varying  size*® 
will  provide  a  more  descriptive  mechanism  for  simulation  of  multiple  (»  3)  condi¬ 
tions. 

8.  SAMPLE  SEQUENCES  OF  TWO  INTERRELATED  WEATHER  CONDITIONS 
8.1  A  48-h  Sequence  of  Ceiling  and  Visibility 

In  temperate  latitudes  and  continental  regions,  the  hour-to-hour  correla¬ 
tion  coefficient,  in  January,  has  been  found  to  be  approximately  prr  =  0.  95  for 

11  ^  ^ 

ceilings,  Pyy  =  0.  92  for  visibilities.  The  correlation  coefficients  between 
ceiling  and  visibility  without  time  lag  have  been  found  to  vary  significantly  with 
the  specific  values  of  ceiling  and  visibility  (Table  3).  This  effect  on  the  simula¬ 
tion  process  is  examined  in  this  section.  Additionally,  since  the  lag  correlation 
coefficients  between  ceiling  and  visibility  should  be  smaller  than  the  non-lag 
correlation  coefficients,  it  is  desirable  to  see  what  happens  when  these  lag 
correlation  coefficients  (p^y,  Py^)  are  *ess  than,  or  greater  than,  the  values 
given  by  the  Whiton  limitation  [Eq.  (16)). 

Figure  7a  shows  the  changes  in  END  of  ceiling  height  for  48  consecutive 
hours.  They  were  generated  stochastically  by  use  of  the  above  equations,  when 
the  intercorrelation  (r^y)  between  ceiling  and  visibility  was  set  at  zero.  That 
is,  ceiling  and  visibility  are  represented  as  changing  independently  of  one 
another.  Figure  7b  shows  the  corresponding  changes  of  visibility.  The  dashed 
lines  on  Figures  7a  and  7b  are  based  on  Bedford,  Mass.  RUSSWO  data  for  the 
month  of  January.  They  indicate  that,  climatically  speaking,  ceiling  is  unlimited 
about  50  percent  of  the  time,  and  visibility  is  unlimited  approximately  60  percent 
of  the  time,  varying  little  with  time  of  day.  As  long  as  the  END  is  greater  than 
that  shown  by  the  upper  dashed  line,  ceiling  or  visibility  is  unlimited.  The 
imaginary  sample  of  Figures  7a  and  7b  is  such  that  the  ceiling  was  virtually 
unlimited  throughout  the  2  days,  but  visibility  was  frequently  restricted  to  less 
than  10  miles. 

Figure  7c  shows  joint  variations  of  ceiling  and  visibility  when  there  is  a 
relatively  small  non-lag  correlation  coefficient  of  0.  4  between  the  END's  of 
ceiling  and  visibility  and  lag  correlation  coefficients  of  0.  34  and  0.  35.  They 
were  produced  by  the  same  random  numbers  as  were  used  in  Figures  7a  and  7b. 

Figures  8(a),  8(b),  and  8(c)  show  a  stochastic  48-h  sequence  of  Bedford, 
Mass. ,  January  ceiling  and  visibility  when  there  was  supposedly  a  non-lag 

10.  Gringorten,  1. 1.  (197  9)  Probability  models  of  weather  conditions  occupying 

a  line  or  an  area,  J.  Appl.  Meteorol.  D!(No.  8):957-977. 

11.  Gringorten,  1. 1.  (1966)  A  stochastic  model  of  the  frequency  and  duration  of 

weather  events,  J.  Appl.  Meteorol.  5(No.  5):606-624. 


Figure  7a.  Simulation  of  END's  of 
Ceiling  Heights  in  a  48-h  Sequence, 
When  There  is  Zero  Correlation 
With  Visibility  (rev  =  0).  Dashed 
lines  show  climatic  frequencies  of 
ceiling  heights  at  Bedford,  Mass, 
in  January.  Hour-to-hour  autocor¬ 
relation  of  ceilings  is  0.  95,  of  vis¬ 
ibility  0.  92 


Figure  7c.  Simulation  of  Ceiling  (Sol¬ 
id  Curve)  and  Visibility  (Broken 
Curve)  When  the  Intercorrelation  is 
rCv  =  0.4  Without  Lag,  pfv  =  0.34, 
PVC  =  °-35  With  1-h  Lag.  Lines  for 
the  climatic  frequencies  are  omitted; 
they  would  be  identical  to  those  in 
Figures  7a  and  7b 


Figure  7b.  Simulation  of  Visibility 
in  a  48-h  Sequence,  When  There  is 
Zero  Correlation  With  Ceiling  (rcv 
=  0).  Dashed  lines  show  climatic  V 
frequencies  of  visibility  at  Bedford, 
Mass,  in  January.  Autocorrela¬ 
tions  are  as  in  Figure  7a 


slightly  less  (p^y  =  0.  34.  -  y(.  =  0.35).  In  Figure  8(e)  .lit-  lag  can  r  lion  coel- 
ficients  were  made  slightly  greater  (Pcv  =  PyC  =  0.35).  Ironically,  the-ceiling 
and  visibility  were  coupled  together  most  closely  when  the  lag  correlation  coef¬ 
ficients  were  the  smallest,  but  this  result  was  not  consistently  repeated  in  other 
trials. 

Figures  9(a),  9(b),  and  9(c)  show  the  48-h  sequence  of  Bedford,  Mass., 
January  ceiling  and  visibility,  stochastically  produced  with  the  same  random 
numbers  as  in  Figures  8(a),  8(b),  and  8(c).  This  time  there  was  supposedly  a 
non-lag  correlation  coefficient  of  0.  95  between  the  END's  of  ceiling  and  visibility. 
The  lag  correlation  coefficients  (p^-y,  Py^.)  were  as  *ow  as  could  be,  for 
practical  solution,  in  Figure  9(a).  The  Whiton  values  were  used  in  Figure  9(b), 
and  relatively  high  values  in  Figure  9(c).  It  appears  that  the  correlation  coef¬ 
ficient  (r^jy)  between  ceiling  and  visibility  made  its  greatest  difference  on  the 
rare  events,  close  to  zero  ceiling  and/or  zero  visibility,  when  the  intercorrela- 
tion  is,  in  fact,  the  highest. 

Faced  with  an  array  of  correlation  coefficients,  such  as  in  Table  3,  the 
easiest  single  value  to  obtain  is  the  average.  In  Table  3,  the  average  is  r^y  = 

0.  86  with  standard  deviation  0.  09.  For 

rCy  =  0.  86 

Pcv  =  (0.  86)(0.  92)  =  0.791 

Pvt  =  (0.  86)(0.  96)  =  0.826 

the  stochastic  generation  of  the  sample  sequence,  with  the  same  random  numbers 
as  in  Figures  8  and  9,  produced  the  sequence  in  Figure  10. 

The  inference  to  be  drawn,  in  the  comparison  of  Figures  8,  9,  and  10,  is 
that  the  average  intercorrelation  r^-y  =  0.86  is  an  acceptable  compromise,  to  be 
used  instead  of  an  assortment  of  values  (Table  3).  If,  however,  the  non-lag 
correlation  coefficient  (r^y)  is  reduced  to  a  single  value,  then  the  lag  correla¬ 
tion  coefficients  (PCy.  Pv(-.)  also  should  be  chosen  for  simplicity,  thus  the 
Whiton  assumption  (Eqs.  16). 

8.2  Simulation  of  a  Two-Station  48-h  Sequence  of  Sky  Covers 

Suppose  that  the  requirement  is  to  simulate  time  changes  of  sky  cover  at  two 
neighboring  stations,  like  Boston  and  Bedford,  Mass,  in  August  when  the  diurnal 
effect  is  large.  Each  autocorrelation,  per  hour  is, 

Pn  =P22  =  0.95  . 

Suppose  the  intercorrelation  is  high,  as  it  probably  is: 
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rl2  =0-95 

P12  ■  (0.95)  (0.  95)  =  0.  9025 

P21  ■  (0.  95)(0.  95)  -  0.9025 
12 

Somerville  and  Bean  have  provided  values  for  the  two  parameters  (a,  0)  in  the 
S -distribution  of  the  cumulative  probability,  F(x),  of  the  sky  cover  x  *  0.0(0. 1) 
1.0.  The  formula  is 

12.  Somerville,  P.N.,  and  Bean,  S.J.  (1979)  A  New  Model  for  Skv  Cover. 
Scientific  Report  No.  5,  Contract  F19628-77-C-0080,  AD  07  8  3  68. 
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Figure  10.  Similar  to  Figure  8  With  Average  Intercorrelations 

F(x)  =  1  -  (1  -  xa)P  .  (37) 

With  the  Somerville  and  Bean  values  used  in  the  program  (Appendix  C),  several 
48-h  sequence  were  obtained.  Figure  11  depicts  a  close  relationship  of  the  two 
sky  covers  that  remain  broken  to  overcast  during  most  of  the  2 -day  period, 
yielding  to  scattered  conditions  during  the  most  likely  time  of  the  evening. 

Figure  12,  with  another  sequence  of  random  numbers,  depicts  clear  to  scattered 
clouds  in  the  first  day,  with  the  likely  development  of  cloud  in  the  afternoon. 
Figure  12  also  depicts  a  10-h  interval  when  the  sky  cover  differed  markedly 
between  stations.  For  the  most  part,  however,  the  changes  are  concurrent. 

Figures  11  and  12  should  be  compared  with  Figure  6.  While  Figures  11  and 
12  give  direct  information  on  sky  cover,  the  plot  of  Figure  6,  on  Normal 
Probability  Paper,  enables  us  to  plot  the  climatic  frequencies  for  comparison 
with  the  synoptic  events.  However,  this  is  only  a  visual  aid.  For  gaming 
purposes,  the  direct  amount  of  sky  cover  may  be  desired  as  the  product  of  the 
computer  exercise. 

9.  DISCUSSION  AND  CONCLUSIONS 

The  simulation  of  a  time  sequence  of  changes  of  two  interrelated  variables 
is  accomplished  by  using  generating  equations,  basically  Eq.  (7).  The  coef¬ 
ficients  (a's  and  b's)  are  determined  in  terms  of  correlation  coefficients  [Eqs. 

( 12 )  and  (13)].  The  terms  symbolized  by  r  are  the  correlation  coefficients 
between  the  variables  without  time  lag,  those  symbolized  by  p  are  correlation 
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Figure  11.  Stochastic  Simulation  of  a  48-h  Sequence  of 
Sky  Cover,  When  Autocorrelations  are  0.  95,  r ^  =  0-  9: 
p  12  =  P21  =  0.9025.  The  curves  are  for  Bedford,  Mass 
and  Boston,  Mass.  August  sky  cover 


Figure  12.  A  Second  Stochastic  Sky-Cover  Simulation. 
Conditions  are  the  same  as  in  Figure  11 


coefficients  with  time  lag.  The  term  p  ,  for  example,  is  the  autocorrelation 
of  ceiling  with  itself  in  time  lag  1  h. 

Unfortunately  the  correlation  coefficients,  lag  correlations,  or  intercorrela 
tions  between  variables,  are  not  readily  available,  especially  between  END's. 
Special  studies,  using  the  hourly  sequences  of  many  years  (10  or  more)  in  each 
month  need  to  be  undertaken,  to  find  such  correlation  coefficients.  A  measure 
of  the  intercorrelation,  without  time  lag,  between  ceiling  and  visibility  is  obtain¬ 
able  from  a  RUSSWO. 

Gringorten*  *  published  hour-to-hour  autocorrelations  of  END's  for  ceilings, 
visibilities  and  other  weather  elements  at  Minneapolis,  Minn.  (Table  4).  There 
were  two  estimates  for  each  correlation;  the  first  estimate  was  best  for  the 
duration  of  the  weather  element  above  a  certain  minimum,  the  second  estimate 
was  best  for  the  duration  of  the  weather  element  below  a  certain  maximum. 
Generally  speaking,  the  hour-to-hour  correlation  coefficient  was  higher  for 
ceiling  than  for  visibility.  Other  studies  have  indicated  higher  autocorrelation 
in  winter  than  in  summer  although  the  figures  for  Minneapolis,  Minn,  do  not 
support  this  conclusion.  A  "best"  a  priori  approximation  for  hour-to-hour  cor¬ 
relation  is  estimated  to  be  0.  95. 


Table  4.  Estimates  of  Hour-to-Hour  Autocorrelations  of  END's  for  Ceil¬ 
ings  and  Visibilities  at  Minneapolis,  Minn. ,  in  Four  Mid-Season  Months 


Element 

Month 

January 

April 

July 

October 

Ceiling 

0.  95 

0.  95 

0.96 

0.  97 

0.  90 

0.  92 

0.  96 

0.  97 

V  isibility 

0.  90 

0.  94 

0.  92 

0.  95 

0.  90 

0.  95 

0.  92 

0.  95 

Table  2  supports  previous  conclusions  that  correlation  coefficients  are  lower 
in  summer  than  in  winter.  Tables  2  and  3  reveal  a  regrettable  dependence  of 
the  tetrachoric  correlation  coefficient  on  the  measure  of  the  ceiling  and  visibility. 
An  overall  average  correlation  coefficient,  however,  works  well,  to  stochastical¬ 
ly  generate  a  sequence  of  values  of  two  interrelated  variables.  Additionally,  a 
time  lapse  simulation  is  well  served  by  the  Whiton  estimates  [Eqs.  (16))  of  lag 
correlations. 
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Appendix  A 

Procedure  to  Generate  a  Simulated  Sequence  of 
END'S  in  an  O-U  Process 


Step  0.  Begin  with  an  arbitrary  initial  END,  [y(t)] .  Choose  a  value  for  the  basic 
time  interval  (6t  hours).  Fix  on  the  parameter:  Relaxation  Time  (  T  hours). 
Decide  on  the  size,  or  the  number  of  END-values  in  sequence  (N). 

Find 


P  =  exp  (-6 1/  T  ) 

Initialize  n  =  0. 

Step  1.  For  n  =  0,  N-l 

Increase  n  by  1.  Use  Subroutine  (A)  to  find  a  random  normal  number, 
n(t  +  6t). 

Find 

y(t  +  6t)  =  p  •  y(t)  +  ^  1  -  p2  .  rjCt  +  6t) 


If  n  <  N,  replace  y(t)  with  y(t  +  <5t)  and  repeat  Step  1. 
If  n  =  N,  the  sequence  is  completed. 


Subroutine  A:  For  generating  a  random  normal  number. 

(Note:  This  is  only  a  simple  suggested  routine.  There  are  others  that  are 
better. ) 

Step  A.  Begin  with  an  arbitrary  random  number  (0  s  xQ  <  l)  or  the  last  number 
in  the  memory. 

Set  J  =  12 

I  x.  =  0. 

Step  B.  For  j  =  0,  J-l 
Increase  j  by  1 

Xj  =  (Xj_j  +  tf)8  -  int^  (Xj.|+  *)8 

Add  x.  to  X  x. 

J  J 

If  j  <  J  replace  with  x^,  and  repeat  Step  B 
For  j  '  J,  find  r)(t  +  fit)  =  £  x.  -  6 
RETURN  to  Step  1. 


^  int  is  the  integral  part  of  a  number,  for  example,  int  (12.375)  =  12 


Appendix  B 

Procedure  to  Generate  a  Simulated  Sequence  of  Two 
Interrelated  ENO't  (y^,  yj)  in  a  Markov  Process 


Step  0.  Assign: 

The  non-lag  correlation  coefficient  between  y1  and  y2 :  r^2 
The  autocorrelation  between  yjft)  and  y^(t  +  6t): 

The  lag  correlation  coefficient  between  y ^ (t)  and  y2(t  +  fit):  p  J2 

The  lag  correlation  coefficient  between  y2  (t)  and  y ^ (t  +  fit):  p2J 

The  autocorrelation  between  y2 <t)  and  y2<t  +  fit):  P22 

Choose  a  value  for  the  time  interval:  fit 

Decide  on  the  number  ofj  pairs  *n  sequence:  N. 


all  =  ^11  "  rl2  p2  1>,/(1  ‘  r12  * 

a  12  =  ^P21  *  rl2  pll>/^1  "  rl2  * 

a21  =  <P12  r  12  p22)/(l  "  1  12  ) 

a22  =  ^P22  "  r  12  p  12  ^  ^ 1  "  1  12  ^ 

hl  =  J1  -  ipn  +P122  *2rl2pnp12)/(1  -  r!22) 


b2  =  \  1  ‘  ^*2  1  +  P22  ”  2rl2P21P22^^  "  rl2  * 

h  =  [r12(l  -  a j  ia22  "  a  12 a2  1  ^  '  ^alla21  +  a12a22^  /blb2 
Begin  with  an  arbitrary  pair  of  initial  random  END's:  y^O),  y2(0). 
Initialize  n  =  0 
Step  1.  For  n  =  0,  N-l 
Increase  n  by  1. 

Use  subroutine  (A)  to  find  a  random  normal  number:  rj^(t  +  fit) 

Use  subroutine  (A)  to  find  a  random  normal  number:  r)2<t  +  fit). 
Find: 

n2(t  +  fit)  =  h  •  r) x <t  +  fit)  +  J 1  -  h2  •  p*(t  +  fit) 


Find: 

yj(t  +  fit)  =  aix  •  y^t)  +  aj2  •  y2 (t)  +  b1  •  Hj(t  +  fit) 

y2(t  +  fit)  =  a21  •  yx<t)  +  a22  •  y2 <t)  +  b2  ■  n2(t  +  fit) 

If  n  <  N,  replace  yj(t)  with  y j (t  +  fit) 

y2(t)  with  y2(t  +  fit)  and  repeat  Step  1. 

If  n  =  N,  the  sequence  is  completed. 

Subroutine  A:  For  generating  a  random  normal  number  (the  same  as  in  Appendix 


Appendix  C 

Procedure  to  Generate  e  Simulated  Sequence  of  Two  Inter¬ 
related  Weather  Element*  (X^,  X2)  in  a  Markov  Procee* 

Step  0.  (Same  as  in  Appendix  B) 

Step  1.  For  n  -  0,  N-l 
Increase  n  by  1. 

Find  a  random  normal  number: 

Find  a  random  normal  number: 

Find: 

r)2(t  +  fit)  =  h  •  rjj  (t  +  fit)  +  ^  1  -  h2  •  rj2(t  +  fit) 

y^t  +  fit)  ■  a^  •  yj(t)  +  al2  •  y2 (t)  +  b.  •  rj^t  +  6t) 

y2(t  +  6t)  =  a21  •  y^t)  +  a22  •  y2(t)  +  bg  •  rj2(t  +  fit)  . 

Step  2.  Use  Subroutine  (B)  to  find  the  probability  (pj)  corresponding  to  yj(t  +  fit) 
Use  Subroutine  (B)  to  find  the  probability  (p2)  corresponding  to  y2(t  +  fit) 


nt(t  +  fit) 
n2(t  +  fit). 


Find  the  hour  of  the  day  (L): 

L  =  24  J  fra  (n/24) }  * 

Find  the  hourly  group  (g): 
g  =  int(L/3)^ 

Use  Subroutine  (C)  to  find  Xj(t  +  fit)  corresponding  to  Pj 
Use  Subroutine  (D)  to  find  (t  +  fit)  corresponding  to  p2 
If  n  <  N,  replace  y^(t)  with  y^(t  +  fit) 

replace  y2(t)  with  y2 (t  +  fit)  and  repeat  Steps  1,2. 

If  n  =  N,  the  sequence  is  completed. 

Subroutine  A:  For  generating  a  random  normal  number  (the  same  as  in  Appendix 
A). 

Subroutine  B:  For  transforming  an  END  (y)  into  probability,  p  =  P(^  y) 

Assign: 

c _  =  0.  196  854 

c2  =  0. 115  194 

c3  =  0.000  344 

c4  =  0.019  527 

Set  l  -  1,  m  =  -1  for  y  2:  0 

Set  I  =  0,  m  =  1  for  y  <  0. 

Find  p  =  l  +  m'|2(l  +  c1ly|+c2y2  +  c3ly3|  +  c4y4)4  | 

RETURN. 

Subroutine  C:  For  ceiling  height  (C),  given  P  (£  C) 

Find,  from  computer  memory,  the  values  a,  b,  c  for  the  hourly  group  (g). 
Where  p  =  P(£  C) 

Find 

C  =  c  [(1  -  p)l/b  -  lj  l/a 
RETURN  or  alternative. 

^fra  is  the  fractional  part  of  a  number,  for  example,  fra(12.  375)  =  0.  37  5 
^int  is  the  integral  part  of  a  number,  for  example,  int(l2.  375)  =  12 


Subroutine  C:  For  cloud  cover  (x),  given  P(^  x) 


Find,  from  computer  memory,  the  values  a,  [3  for  the  hourly  group  (g) 
Where  p  =  P(s  x) 

Find 

x  =  [  1  -  (1  -  p)1^]  1^a 
RETURN. 

Subroutine  D:  For  visibility  (V)  given  PCs  V) 

Find,  from  computer  memory,  the  values  a,  ji  for  the  hourly  group  (g) 
Where  p  =  P(s=  V) 

Find 


RETURN,  or  alternative. 


Subroutine  D:  For  cloud  cover  (same  as  for  Subroutine  C  for  cloud  cover). 
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